Reliability-based topology optimization design method for part structure considering bounded hybrid uncertainties

ABSTRACT

A reliability-based topology optimization design method for a part structure considering bounded hybrid uncertainties, includes the following steps: considering the uncertainties in the manufacture and service of the part structure, describing an external load with insufficient samples and a material property with sufficient samples as an interval variable and a bounded probabilistic variable respectively; discretizing a design domain of the part structure, setting the physical and geometric constraints, and establishing a reliability-based topology optimization design model; solving by a moving asymptote algorithm: decoupling the probabilistic and interval uncertainties, and determining the worst working condition by using gradients of constraint performance functions; defining a performance fluctuation under the worst working condition and calculating reliability of constraint performance; and finally, calculating the gradients of objective and constraint functions with respect to the design variables for iteration.

CROSS-REFERENCE TO RELATED APPLICATIONS

The present application is a continuation of International Application No. PCT/CN2021/079563, filed on Mar. 8, 2021, the content of which is incorporated herein by reference in its entirety.

TECHNICAL FIELD

The present disclosure belongs to the field of the optimization design of equipment structure, and relates to a reliability-based topology optimization design method for a part structure considering bounded hybrid uncertainties.

BACKGROUND

As a method of distributing finite materials in a design domain to optimize the structural objective performance, topology optimization has been widely applied in product design. Since there are various uncertainties in the processes of manufacture and use, in order to ensure that theoretical results obtained by topology optimization will not deteriorate in performance after actual manufacture, the impact of the uncertainties must be considered in the design stage. However, traditional reliability-based topology optimization method that uses only a single type of mathematical model to describe multi-source uncertainties may probably lead to the distortion of uncertainty modeling and the failure of the resulting optimal scheme of a product structure, so it is necessary to carry out the reliability-based topology optimization of a part structure simultaneously considering probabilistic and interval hybrid uncertainties.

In addition, in the field of reliability-based structural topology optimization considering hybrid uncertainties at present, most of researches still transform an original problem into an equivalent deterministic problem after searching the most probable point, which has obvious defects for an optimization problem involving multiple reliability-based constraints: the most probable points obtained from the analyses of different reliability-based constraints are often inconsistent, and their trade-offs either introduce subjective factors or involve additional calculations. In addition, the reliability expressions constructed by existing methods often requires additional condition discrimination, and the calculation result of reliability-based gradient in some scenarios can hardly apply to a gradient-based topology optimization framework. Therefore, on the basis of giving the description of bounded hybrid uncertainties, it is necessary to study a reliability-based topology optimization method involving multiple constraint performance functions, and improve the gradient analyticity of the reliability expressions of constraint performance.

SUMMARY

In order to solve the problem of reliability-based topology optimization design of a part structure under the influence of multi-source uncertainties, the present disclosure provides a reliability-based topology optimization design method for a part structure by considering bounded hybrid uncertainties. The method includes the following steps: considering the uncertainties in the manufacture and use of a part structure, describing an external load with insufficient samples as an interval uncertainty and a material property with sufficient samples as a bounded probabilistic uncertainty; discretizing a design domain, setting physical and geometric constraints, and establishing a reliability-based topology optimization design model of the part structure; iteratively solving by a moving asymptote algorithm: decoupling the probabilistic and interval uncertainties, and determining a worst working condition by the gradients of constraint performance functions; defining a performance fluctuation under the worst working condition and calculating the reliability of constraint performance; and finally, calculating the gradients of objective and constraint functions on design variables for iteration. The present disclosure efficiently solves the problem of the reliability-based topology optimization design of a part structure when probabilistic and interval uncertainty factors coexist, and is very valuable to engineering application.

The present disclosure is implemented through the following technical solution: a reliability-based topology optimization design method for a part structure by considering bounded hybrid uncertainties includes the following steps:

-   -   1) Considering following uncertainties in the manufacture and         service processes of a part structure: regarding the amplitude         and direction of an external load which are difficult to obtain         sufficient sample information as interval uncertainties; and         regarding a material property of the part structure with         sufficient sample information as a bounded probabilistic         uncertainty, and describing a bounded probabilistic uncertainty         parameter as a random variable subject to generalized beta         distribution.     -   2) Discretizing a design domain of the part structure, which is         specifically as follows:

Simplifying a force condition of the part structure into a two-dimensional plane stress state, retaining installation holes and removing structural details to improve calculating efficiency; placing the simplified part structure in a regular rectangular design domain, and dividing the rectangular design domain into N_(x)×N_(y) square elements, where N_(x) and N_(y) are the numbers of divisions along x,y axes respectively; and based on a solid isotropic material with penalization (SIMP) topology optimization framework, imparting each element with a unique design variable ρ_(e)∈[0,1] (e=1, 2, . . . , N_(x)·N_(y)).

-   -   3) Imposing physical constraints and geometric constraints on         the discretized structure, which is specifically as follows:     -   3.1) Imposing the physical constraints including fixing or         supporting and an external load according to a classical finite         element method.     -   3.2) Imposing the geometric constraints including specified         holes in the structure and areas where materials are forcibly         retained.

A method is to set design variables corresponding to elements in the holes as ρ_(e)≡0 and set design variables corresponding to elements in the area where the materials are to be retained as ρ_(e)≡1, and values thereof are kept unchanged in subsequent optimization process.

-   -   4) By taking a space utilization rate of the design domain as an         objective function and taking the displacements of several key         points of the part structure considering a joint influence of         interval and bounded probabilistic hybrid uncertainties as         reliability constraint performance, establishing a         reliability-based topology optimization design model for the         part structure, as shown in Eq.1:

$\begin{matrix} {{\min\limits_{\rho}{V(\rho)}} = {\sum\limits_{e = 1}^{N_{e}}{\rho_{e}/V_{0}}}} & {{Eq}\text{.1}} \end{matrix}$ s.t.g_(q)(ρ, X, I) = −P(u_(q)(ρ, X, I) ≤ u_(qcri)) + P_(q) ≤ 0(q = 1, 2, ⋯, N_(con)) K(ρ, X)U = F(I) 0 ≤ ρ_(min) ≤ ρ_(e) ≤ 1(e = 1, 2, ⋯, N_(e)) X = (X₁, X₂, ⋯, X_(m))^(T), I = (f₁, f₂, ⋯, f_(n), α₁, α₂, ⋯, α_(n))^(T)

-   -    where in Eq.1, ρ=(ρ₁, ρ₂, . . . , ρ_(N) _(x) _(·N) _(y) )^(T)         is a design vector composed of the design variables ρ_(e)(e=1,         2, . . . , N_(x)·N_(y)), and ρ_(min) is a minimum allowable         value for the design variables, and there is ρ_(min)=1E−6; the         total number of elements is N_(e)=N_(x)·N_(y); a bounded         probabilistic uncertainty vector X=(X₁, X₂, . . . , X_(m))^(T)         comprises m uncertain material properties of the part structure;         an interval uncertainty vector I=(f₁, f₂, . . . , f_(n), α₁, α₂,         . . . , α_(n))^(T) comprises amplitudes f₁, f₂, . . . , f_(n)         and direction angles α₁, α₂, . . . , α_(n) of n uncertain         external loads on the part structure.

V(ρ) is the space utilization rate of the design domain, corresponding to a total material usage of the part structure, and V₀ is a volume of the design domain.

g_(q)(ρ,X,I) is the q^(th) constraint function; u_(q)(ρ,X,I) is the displacement at the q^(th) key point regarded as the q^(th) constraint performance, which is a function of the design variables and the uncertainty variables, and for simplicity in the following, it is briefly recorded as u_(q), and u_(qcri) is a maximum allowable value thereof; P(⋅) calculates an occurrence probability of an event in a bracket, P_(q) is a reliability index of the q^(th) constraint performance, and N_(con) is the number of constraint functions.

In an equilibrium equation K(ρ,X)U=F(I) of the part structure, K(ρ,X) is a (2(N_(x)+1)(N_(y)+1))×(2(N_(x)+1)(N_(y)+1))-dimensional global stiffness matrix, which is affected by the design vector p and the bounded probabilistic uncertainty vector X; F(I) is a (2(N_(x)+1)(N_(y)+1))-dimensional nodal force vector, which is affected by the interval uncertainty vector I; U is a (2(N_(x)+1)(N_(y)+1))-dimensional nodal displacement vector; and u_(q) is extracted from U according to u_(q)=L_(q) ^(T)U.

In the (2(N_(x)+1)(N_(y)+1))-dimensional column vector L_(q), except that the element at the position corresponding to the q^(th) key point is 1, other elements are all 0.

-   -   5) Calculating the reliability of a constraint performance under         bounded hybrid uncertainties:     -   5.1) Searching the worst working condition of constraint         performance u_(q) first:     -   5.1.1) Setting X=μ_(X)=(μ_(X) ₁ , μ_(X) ₂ , . . . , μ_(X) _(m)         ), where μ_(X) ₁ , μ_(X) ₂ , . . . , μ_(X) _(m) are the means of         the uncertainties X₁, X₂, . . . , X_(m); under the circumstance,         the constraint performance u_(q) is only affected by the         interval uncertainty; and an uncertain external load F_(s) (s=1,         2, . . . , n) is rewritten as its components in the horizontal         and vertical directions such as F_(s)=[f_(s) cos α_(s), f_(s)         sin α_(s)]^(T).

5.1.2) Based on a linear elastic hypothesis in small deformation, displacement U_(s) caused by F_(s) is calculated by Eq.2 based on e_(x) ^(x)=[1, 0]^(T), e_(s) ^(y)=[0 1]^(T) and F_(s), where e_(s) ^(x) and e_(s) ^(y) are respectively the element nodal forces in the horizontal and vertical directions at the point where F_(s) exerts;

U _(s) =U _(s) ^(x) +U _(s) ^(y) =u _(s) ^(x) f _(s) cos α_(s) +u _(s) ^(y) f _(s) sin α_(s) =[u _(s) ^(x) ,u _(s) ^(y) ]·F _(s)  Eq.2

In Eq. 2, u_(s) ^(x)=[u_(s) ^(x) 0]^(T) and u_(s) ^(y)=[0 u_(s) ^(y)]^(T) are nodal displacement vectors calculated through the equilibrium equation of the part structure when only the elemental nodal force e_(s) ^(x) or e_(s) ^(y) acts, respectively.

5.1.3) Calculating the gradients of the constraint performance u_(q) with respect to the amplitude and direction of an uncertain external load according to Eq.3 and Eq.4:

$\begin{matrix} {\frac{\partial u_{q}}{\partial f_{s}} = {{L_{q}^{T}\frac{\partial U}{\partial f_{s}}} = {{L_{q}^{T}{\sum\limits_{i = 1}^{n}\frac{\partial U_{i}}{\partial f_{s}}}} = {L_{q}^{T} \cdot \left\lbrack {u_{s}^{x},u_{s}^{y}} \right\rbrack \cdot \left\lbrack {{\cos\alpha_{s}},{\sin\alpha_{s}}} \right\rbrack^{T}}}}} & {{Eq}\text{.3}} \end{matrix}$ $\begin{matrix} {\frac{\partial u_{q}}{\partial\alpha_{s}} = {{L_{q}^{T}\frac{\partial U}{\partial\alpha_{s}}} = {{L_{q}^{T}{\sum\limits_{i = 1}^{n}\frac{\partial U_{i}}{\partial\alpha_{s}}}} = {L_{q}^{T} \cdot \left\lbrack {u_{s}^{x},u_{s}^{y}} \right\rbrack \cdot \left\lbrack {{{- f_{s}}\sin\alpha_{s}},{f_{s}\cos\alpha_{s}}} \right\rbrack^{T}}}}} & {{Eq}\text{.4}} \end{matrix}$

-   -   5.1.4) By utilizing the results of Eq.3 and Eq.4, solving the         worst working condition Ĩ_(q) by a gradient search algorithm,         and the external load of the worst working condition is F={tilde         over (F)}_(q) at the moment.     -   5.2) Restoring μ_(X) to X under the worst working condition         Ĩ_(q), so the constraint performance ũ_(q)=u_(q) (ρ,X,Ĩ_(q))         under the worst working condition is only manifested as a         probabilistic type, and solving the performance fluctuation         under the worst working condition to evaluate reliability,         details are as follows:     -   5.2.1) Calculating, according to Eq.5, a gradient of the         constraint performance ũ_(q) with respect to a bounded         probabilistic uncertainty parameter X_(i) (i=1, 2, . . . , m)         under the worst working condition:

$\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{q}}{\partial X_{i}} = {{L_{q}^{T}\frac{\partial\overset{\sim}{U}}{\partial X_{i}}} = {{{- L_{q}^{T}}K^{- 1}\frac{\partial K}{\partial X_{i}}\overset{\sim}{U}} = {{- L_{q}^{T}}{K^{- 1}\left( {\sum\limits_{e = 1}^{N_{e}}{\rho_{e}^{p}\frac{\partial k_{e}}{\partial X_{i}}}} \right)}\overset{\sim}{U}}}}} & {{Eq}\text{.5}} \end{matrix}$

-   -    where a summation symbol is a combination operation of the         element stiffness matrix defined by a finite element theory,         k_(e) is an element stiffness matrix, the nodal displacement         vector Ũ under the worst working condition is obtained by         solving the governing equation KŨ={tilde over (F)}_(q), and the         penalty factor is generally p=3.     -   5.2.2) Based on the result of Eq.5, searching a bounded         probabilistic uncertainty vector minimizing or maximizing u_(q)         (ρ,X,Ĩ_(q)) respectively according to Eq.6:

$\begin{matrix} \left\{ {\begin{matrix} {X_{q}^{L} = {\underset{X}{\arg\min}{u_{q}\left( {\rho,X,{\overset{\sim}{I}}_{q}} \right)}}} \\ {X_{q}^{R} = {\underset{X}{\arg\max}{u_{q}\left( {\rho,X,{\overset{\sim}{I}}_{q}} \right)}}} \end{matrix}\left( {{q = 1},2,\cdots,N_{con}} \right)} \right. & {{Eq}\text{.6}} \end{matrix}$

-   -    where corresponding to X_(q) ^(L) and X_(q) ^(R) the global         stiffness matrix K is K_(q) ^(L) and K_(q) ^(R) respectively,         which are uniformly marked as K_(q) ^(*) (*=L,R), and the nodal         displacement vector Ũ is Ũ_(q) ^(L) and Ũ_(q) ^(R) respectively,         which are uniformly marked as Ũ_(q) ^(*) (*=L,R); and     -   5.2.3) Briefly recording ũ_(q) ^(L)=u_(q)(ρ,X_(q) ^(L),Ĩ_(q))         and ũ_(q) ^(R)=u_(q)(ρ,X_(q) ^(R),Ĩ_(q)), defining [ũ_(q)         ^(L),ũ_(q) ^(R)] as the fluctuation of the constraint         performance u_(q) under the worst working condition, and         calculating the reliability {tilde over (R)}_(q) of the         constraint performance according to Eq.7:

$\begin{matrix} {{\overset{\sim}{R}}_{q} = {{\frac{1}{2}\tan h\left\{ {P \cdot \left( {u_{qcri} - u_{q}^{C}} \right) \cdot \left\lbrack {1 + \left( {P \cdot \left( {u_{qcri} - u_{q}^{C}} \right)} \right)^{\gamma}} \right\rbrack} \right\}} + \frac{1}{2}}} & {{Eq}\text{.7}} \end{matrix}$

-   -    where u_(q) ^(C)=(ũ_(q) ^(R)+ũ_(q) ^(L))/2 is a midpoint of the         performance fluctuation under the worst working condition; a         multiplier is P=1/(u_(q) ^(W)−ε^(u)), where ε^(u) is a small         constant for adjusting the reliability at the boundary position         of the performance fluctuation under the worst working         condition, the value of which is recommended as         ε^(u)=0.00˜0.05u_(q) ^(W) based on actual tests; u_(q)         ^(W)=(ũ_(q) ^(R)−ũ_(q) ^(L))/2 is a radius of the performance         fluctuation under the worst working condition; γ∈{2i|i∈N⁺} is a         regulatory factor.

The constraint function in Eq.1 can be rewritten as Eq.8:

g _(q)(ρ,X,I)=−{tilde over (R)} _(q) +P _(q)≤0 (q=1,2, . . . ,N _(con))  Eq.8

-   -   6) Calculating the gradients of the objective and constraint         functions with respect to the design variables:     -   6.1) Calculating the gradient of the objective function through         Eq.9:

$\begin{matrix} {\frac{\partial{V(\rho)}}{\partial\rho_{e}} = {\frac{1}{V_{0}}\left( {{e = 1},2,\cdots,N_{e}} \right)}} & {{Eq}\text{.9}} \end{matrix}$

-   -   6.2) Solving the gradients of the constraint functions as         follows:     -   6.2.1) Writing a gradient expression of g_(q)(ρ,X,I) according         to a chain rule, as shown in Eq.10:

$\begin{matrix} {\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial\rho_{e}} = {{\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{q}^{L}}\frac{\partial{\overset{\sim}{u}}_{q}^{L}}{\partial\rho_{e}}} + {\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{q}^{R}}\frac{\partial{\overset{\sim}{u}}_{q}^{R}}{\partial\rho_{e}}\left( {{e = 1},2,\cdots,N_{e}} \right)}}} & {{Eq}\text{.10}} \end{matrix}$

-   -   6.2.2) Describing the function in the bracket of tanh(⋅) in Eq.7         as R (u_(qcri)), the gradient terms ∂g_(q)(ρ,X,I)/∂ũ_(q) ^(L)         and ∂g_(q)(ρ,X,I)/∂ũ_(q) ^(R) in Eq.10 can be calculated         according to Eq.11:

$\begin{matrix} {\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{q}^{*}} = {{\frac{{sech}^{2}\left( {R\left( u_{qcri} \right)} \right)}{2} \cdot \frac{\partial{R\left( u_{qcri} \right)}}{\partial{\overset{\sim}{u}}_{q}^{*}}}{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}\text{.11}} \end{matrix}$

-   -    where the gradient term ∂R(u_(qcri))/∂ũ_(q) ^(*) is         specifically as follows:

$\begin{matrix} {\frac{\partial{R\left( u_{qcri} \right)}}{\partial{\overset{\sim}{u}}_{q}^{L}} = {\frac{u_{qcri} - u_{q}^{C}}{2P^{2}} - \frac{P}{2} + \frac{\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma + 1}}{2P^{2}} - {\frac{P\left( {\gamma + 1} \right)}{2}\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma}}}} & {{Eq}\text{.12}} \end{matrix}$ $\begin{matrix} {\frac{\partial{R\left( u_{qcri} \right)}}{\partial{\overset{\sim}{u}}_{q}^{R}} = {{- \frac{u_{qcri} - u_{q}^{C}}{2P^{2}}} - \frac{P}{2} - \frac{\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma + 1}}{2P^{2}} - {\frac{P\left( {\gamma + 1} \right)}{2}\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma}}}} & {{Eq}\text{.13}} \end{matrix}$

-   -   6.2.3) Giving the gradient terms ∂ũ_(q) ^(L)/∂ρ_(e) and ∂ũ_(q)         ^(R)/∂ρ_(e) in Eq.10 in a uniform form according to the SIMP         framework:

$\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{q}^{*}}{\partial\rho_{e}} = {{- {L_{q}^{T}\left( K_{q}^{*} \right)}^{- 1}}\frac{\partial K_{q}^{*}}{\partial\rho_{e}}{\overset{\sim}{U}}_{q}^{*}{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}\text{.14}} \end{matrix}$

-   -    where K_(q) ^(*) and Ũ_(q) ^(*) (*=L, R) are defined in 5.2.2;         and the gradient term ∂K_(q) ^(*)/∂ρ_(e) is calculated according         to Eq.15:

$\begin{matrix} {\frac{\partial K_{q}^{*}}{\partial\rho_{e}} = {p\rho_{e}^{p - 1}\left\langle k_{eq}^{*} \right\rangle{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}\text{.15}} \end{matrix}$

-   -    where k_(eq) ^(*) is the element stiffness matrix extracted         from K_(q) ^(*), and         k_(eq) ^(*)         is a square matrix reconstructed by performing a combined         operation on the elements in k_(eq) ^(*) according to the         element stiffness matrix, and is consistent with K_(q) ^(*) in         dimensionality.     -   6.2.4) Substituting all the gradient terms in Eq.11 to Eq.15         into Eq.10 to obtain the gradient of the constraint function         g_(q)(ρ,X,I).     -   7) Updating the design variables by using a moving asymptote         algorithm based on the gradients of the objective and constraint         functions with regard to the design variables; and checking the         difference between the objective function values in the current         and previous iterations, for the first iteration, defining the         difference as the objective function value, and if the         difference is less than a convergence threshold, outputting the         updated design variables, otherwise repeating steps 5) to 7).

The present disclosure has the beneficial effects:

-   -   1) Following uncertainties in the manufacture and service         processes of the part structure are considered: material         properties of the part structure, and the amplitude and         direction of the external load. Since sufficient sample         information of the external load is difficult to obtain, the         uncertainties of amplitude and direction are regarded as the         interval uncertainties; and the property of a base material with         sufficient sample information is regarded as the bounded         probabilistic uncertainty, which is described as a random         variable subject to generalized beta distribution, so that the         defect of only considering probabilistic or interval uncertainty         in existing reliability-based topology optimization design         method is overcome, and the constructed reliability-based         topology optimization design model for the part structure better         conforms to engineering practice.     -   2) An explicit expression of the constraint performance with         respect to design variables and uncertainty parameters is         established by a classic finite element framework; and the         hypothesis of linear elastic deformation is introduced, the         final deformation of the structure is obtained by superposing         the deformation produced by the separate action of every         external load, and then the gradient information of the         constraint performance with regard to an uncertain external load         is calculated, so that the worst working condition corresponding         to the worst constraint performance of the structure is         obtained, and a theoretical basis is provided for guaranteeing         the safe service of the structure.     -   3) The performance fluctuation under the worst working condition         is defined, and a mathematical expression for calculating the         reliability of the constraint performance under the hybrid         uncertainty condition is given; existing reliability-based         topology optimization methods considering multiple constraint         performance have the defects that the most probable points for         the failures of multiple constraint performance are subjectively         chosen, an uncertainty covariance-correlation matrix needs to be         solved additionally, and reliability expressions usually contain         conditional judgment; and compared with the above existing         methods, the method of the present disclosure fundamentally         avoids solving the most probable points, and investigates the         complete ranges of constraint performance by using the         boundedness of the hybrid uncertainties, no additional         conditional judgment is required by the mathematical expression         for the reliability of constraint performance, and thus the         efficient solution of the gradient information with regard to         design variables is guaranteed.

BRIEF DESCRIPTION OF DRAWINGS

FIG. 1 is a flow diagram of reliability-based topology optimization design for a part structure considering bounded hybrid uncertainties.

FIG. 2 is a three-dimensional external view of a certain type of tunnel boring mechanism and a schematic diagram of a support structure in the outer cutterhead.

FIG. 3 is an initial design diagram for a support structure in the outer cutterhead.

FIG. 4 is a schematic diagram of the design domain for the reliability-based topology optimization of a support structure in the outer cutterhead.

FIG. 5 is the result of the reliability-based topology optimization design of a support structure in the outer cutterhead.

FIG. 6 is a final design diagram of a support structure in the outer cutterhead obtained by smoothing the result of the reliability-based topology optimization design.

DESCRIPTION OF EMBODIMENTS

The present disclosure is further described below in conjunction with the accompanying drawings and specific example.

The information involved in the drawings is the actual application data in reliability-based topology optimization design for a support structure in the outer cutterhead of a certain type of tunnel boring mechanism in the present disclosure, and FIG. 1 is a flow diagram of reliability-based topology optimization design for a part structure considering bounded hybrid uncertainties.

-   -   1. By taking the outer cutterhead support structure in a certain         type of tunnel boring mechanism which is made of a high-strength         low-alloy steel material as shown in FIG. 2 as the research         object, the uncertainties of the support structure in the         processes of manufacture and service are considered:     -   1.1) FIG. 3 illustrates the relevant dimensions of an initial         design of the outer cutterhead support structure of the         tunneling boring mechanism, and FIG. 4 is a boundary condition         set for the reliability-based topology optimization. The support         structure bears an axial load at its top transmitted by the         outer cutterhead in the cutting movement process of the         tunneling boring mechanism; since the support structure is flat,         the optimization problem is simplified into a problem in a         two-dimensional plane, and the load is considered as one         uniformly distributed linear load herein, with its amplitude and         direction being uncertain due to the fluctuation of physical         properties of the rock strata being cut; but it is difficult to         measure the external load in the working process of the         tunneling boring mechanism, and it is difficult to obtain         sufficient sample information of the external load, thus its         amplitude f and direction angle α are regarded as interval         uncertainties.     -   1.2) The Young's modulus E_(M) and Poisson ratio v_(M) of         high-strength low-alloy steel used by the support structure in         the outer cutterhead have relatively obvious uncertainties due         to the non-uniformity of raw material physical properties, the         fluctuation of a metallurgical process and so on, but sufficient         sample information can be obtained by measuring finished         products, so they can be regarded as bounded probabilistic         uncertainties and described as random variables subject to         generalized beta distributions; and the uncertainty information         is summarized in Table 1.

TABLE 1 Summary of the uncertainty information of the support structure in the outer cutterhead of the tunneling boring mechanism Uncertainty Type of uncertainty variable Value range Uncertain parameter* E_(M) (GPa) Bounded [200.00, 210.00] μ_(EM) = 206.00, probabilistic variable σ_(EM) = 1.20 α_(EM) = 5.30, β_(EM) = 6.28 v_(M) Bounded [0.28, 0.32] μ_(vM) = 0.30, probabilistic variable σ_(vM) = 5.00E−3 α_(vM) = β_(vM) = 5.32 f (kN/m) Interval variable [85.4, 87.4] <86.4, 1.0> A Interval variable [−105.00°, −75.00°]  <−90.00°, 15.00°> *Midpoint and radius for interval variables; mean value and standard deviation for bounded probabilistic variables.

-   -   2. The design domain of the support structure is discretized,         which is specifically as follows:

The support structure in the outer cutterhead of the tunneling boring mechanism is flat, so its force condition is simplified into a two-dimensional plane stress state; the support structure to be optimized is placed in a regular rectangular design domain (a range defined by an outermost solid line in FIG. 4 , with the dimension being X×Y=1900 mm×1400 mm), the rectangular design domain is divided into N_(x)×N_(y) square elements, where N_(x) and N_(y) are the numbers of divisions along x, y axes, respectively, which are N_(x)=190 and N_(y)=140 in the design; and based on the classic SIMP framework for topology optimization, each element is imparted with a unique design variable ρ∈[0,1] (e=1, 2, . . . , 190×140)

-   -   3. According to the classic finite element method, physical and         geometric constraints are imposed on the discretized support         structure of the outer cutterhead, which is specifically as         follows:     -   3.1) Physical constraints: all elements at the bottom of the         support structure in FIG. 4 are set as fixed constraints, and         displacement in the y direction is allowed at both left and         right sides; and a uniformly distributed linear load is applied         at the top of the support structure in FIG. 4 , which has an         uncertain amplitude f and an uncertain direction angle α.     -   3.2) Geometric constraints: as shown in FIG. 4 , slash regions         in the design domain are non-design areas, setting design         variables corresponding to the elements in the non-design areas         as ρ_(e)≡0, and keeping their values unchanged in subsequent         optimization process.     -   4. A space utilization rate of the design domain is taken as an         objective function and it is considered that displacements         u_(xA) and u_(yA) in x and y directions of point A (552, 413)         which deforms remarkably in both a normal working condition and         a critical state of the outer cutterhead support structure under         the joint influence of interval and bounded probabilistic hybrid         uncertainties are reliability constraint performance; allowable         values are u_(xAcri)=0.4 mm and u_(yAcri)=1.0 mm, respectively;         and a reliability-based topology optimization design model of         the support structure in the outer cutterhead is established         under the reliability requirement P_(xA)=P_(yA)=0.98, which is         as shown in Eq.16:

$\begin{matrix} {{\underset{\rho}{\min}{V(\rho)}} = {\sum\limits_{e = 1}^{190 \times 140}{\rho_{e}/\left( {190 \times 140} \right)}}} & {{Eq}\text{.16}} \end{matrix}$ s.t.g_(xA)(ρ, X, I) = −P(u_(xA)(ρ, X, I) ≤ u_(xAcri)) + P_(xA) ≤ 0 g_(yA)(ρ, X, I) = −P(u_(yA)(ρ, X, I) ≤ u_(yAcri)) + P_(yA) ≤ 0 K(ρ, X)U = F(I)0 ≤ ρ_(min) ≤ ρ_(e) ≤ 1(e = 1, 2, ⋯, 190 × 140) X = (E_(M), v_(M))^(T), I = (f, α)^(T)

-   -    where ρ=(ρ₁, ρ₂, . . . , ρ_(190×140))^(T) is a design vector,         and the minimum allowable value of all design variables is         ρ_(min)=0.001; X=(E_(M),v_(M))^(T) is a bounded probabilistic         uncertainty vector; I=(f,α)^(T) is an interval uncertainty         vector; V(ρ) is the space utilization rate in the design domain         corresponding to the total material usage; g_(xA)(ρ,X,I) and         g_(yA)(ρ,X,I) are constraint functions; constraint performance         u_(xA)(ρ,X,I) and u_(yA)(ρ,X,I) are the functions of the design         variables and the uncertainty variables, which are briefly         recorded as u_(xA) and u_(yA) in the following for the purpose         of simplicity.

K(ρ,X)U=F(I) K(ρ,X) is an equilibrium equation, where K(ρ,X) is a 2(191×141)×2(191×141)-dimensional global stiffness matrix; F(I) is a 2(191×141)-dimensional nodal force vector; and U is a 2(191×141)-dimensional nodal displacement vector.

-   -   5. Reliability of the constraint performance under the hybrid         uncertainties are calculated, which is explained in detail by         taking g_(xA)(ρ,X,I) as an example in the following:     -   5.1) The worst working condition of constraint performance         u_(xA) is searched first:     -   5.1.1) Setting X=μ_(X)=(μ_(E) _(M) ,μ_(v) _(M) )^(T), where         μ_(E) _(M) ,μ_(v) _(M) are the means of the uncertainties         E_(M),v_(M). Under the circumstance, the constraint performance         u_(xA) is only affected by the interval uncertainty; and an         uncertain external load is rewritten as its components in         horizontal and vertical directions such as F=[f cos α, f sin         α]^(T).     -   5.1.2) Based on a linear elastic hypothesis in small         deformation, displacement U caused by F is calculated by Eq.17         based on e^(x)=[1,0]^(T), e^(y)[0,1]^(T) and F, where e^(x) and         e^(y) are respectively the element nodal forces in the         horizontal and vertical directions at the point where F exerted;

U=U ^(x) +U ^(y) =u ^(x) f cos α+u ^(y) f sin α=[u ^(x) ,u ^(y) ]·F  Eq.17

-   -    where u_(x)=[u^(x),0]^(T) and u^(y)=[0,u^(y)]^(T) are         respectively the nodal displacement vectors calculated through         the equilibrium equation when only the element nodal force e^(x)         or e^(y) acts;     -   5.1.3) The gradients of the constraint performance with respect         to the amplitude and direction of an uncertain external load are         calculated according to Eq.18 and Eq.19:

$\begin{matrix} {\frac{\partial u_{xA}}{\partial f} = {L_{xA}^{T} \cdot \left\lbrack {u^{x},u^{y}} \right\rbrack \cdot \left\lbrack {{\cos\alpha},{\sin\alpha}} \right\rbrack^{T}}} & {{Eq}\text{.18}} \end{matrix}$ $\begin{matrix} {\frac{\partial u_{xA}}{\partial\alpha} = {L_{xA}^{T} \cdot \left\lbrack {u^{x},u^{y}} \right\rbrack \cdot \left\lbrack {{{- f}\sin\alpha},{f\cos\alpha}} \right\rbrack^{T}}} & {{Eq}\text{.19}} \end{matrix}$

-   -   5.1.4) By utilizing the results of Eq.18 and Eq.19, the worst         working condition Ĩ_(xA) is solved by the gradient search         algorithm, and the external load of the worst working condition         is F={tilde over (F)}_(xA) at the moment.     -   5.2) μ_(X) is restored to X under the worst working condition         Ĩ_(xA), so the constraint performance ũ_(xA)=u_(xA)(ρ,X,Ĩ_(xA))         under the worst working condition is only manifested as a         probabilistic type, and the performance fluctuation under the         worst working condition is solved to evaluate reliability,         details being as follows:     -   5.2.1) According to Eq.20 and Eq.21, the gradients of the         constraint performance ũ_(xA) with respect to the bounded         probabilistic uncertainties E_(M),v_(M) under the worst working         condition Ĩ_(xA) are calculated:

$\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{xA}}{\partial E_{M}} = {{- L_{xA}^{T}}{K^{- 1}\left( {\sum\limits_{e = 1}^{190 \times 140}{\rho_{e}^{3}\frac{\partial k_{e}}{\partial E_{M}}}} \right)}\overset{\sim}{U}}} & {{Eq}\text{.20}} \end{matrix}$ $\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{xA}}{\partial v_{M}} = {{- L_{xA}^{T}}{K^{- 1}\left( {\sum\limits_{e = 1}^{190 \times 140}{\rho_{e}^{3}\frac{\partial k_{e}}{\partial v_{M}}}} \right)}\overset{\sim}{U}}} & {{Eq}\text{.21}} \end{matrix}$

-   -    where a summation symbol is a combined operation of element         stiffness matrix defined by the finite element theory, k_(e) is         an element stiffness matrix, and the nodal displacement vector Ũ         under the worst working condition is obtained by solving the         governing equation KŨ={tilde over (F)}_(xA); and the penalty         factor is p=3;     -   5.2.2) Based on results of Eq.20 and Eq.21, searching a bounded         probabilistic uncertainty vector minimizing or maximizing         ũ_(xA)=u_(xA)(ρ,X,Ĩ_(xA)) respectively according to Eq.22:

$\begin{matrix} \left\{ \begin{matrix} {X_{xA}^{L} = {{\underset{X}{\arg\min}{u_{xA}\left( {\rho,X,{\overset{\sim}{I}}_{xA}} \right)}} = \left\lbrack {200.05,0.287} \right\rbrack^{T}}} \\ {X_{xA}^{R} = {{\underset{X}{\arg\max}{u_{xA}\left( {\rho,X,{\overset{\sim}{I}}_{xA}} \right)}} = \left\lbrack {210.,0.309} \right\rbrack^{T}}} \end{matrix} \right. & {{Eq}\text{.22}} \end{matrix}$

-   -    where corresponding to X_(xA) ^(L) and X_(xA) ^(R), the global         stiffness matrix K is K_(xA) ^(L) and K_(xA) ^(R) respectively,         which are uniformly marked as K_(xA) ^(*)(* L, R), and the nodal         displacement vector Ũ is Ũ_(xA) ^(L) and Ũ_(xA) ^(R)         respectively, which are uniformly marked as Ũ_(xA) ^(*)(*=L,R);         and     -   5.2.3) ũ_(xA) ^(L)=u_(xA)(ρ,X_(xA) ^(L),Ĩ_(xA))=0.492 mm and         ũ_(xA) ^(R)=u_(xA)(ρ,X_(xA) ^(R),Ĩ_(xA))=0.533 mm are briefly         recorded, [ũ_(xA) ^(L),ũ_(xA) ^(R)]=[0.492,0.533] mm is defined         as the fluctuation of the constraint performance u_(xA) under         the worst working condition; ε^(u)=0 and γ=4 are set, and the         reliability {tilde over (R)}_(xA) of the constraint performance         is calculated according to Eq.23:

$\begin{matrix} {{\overset{\sim}{R}}_{xA} = {{{\frac{1}{2}\tan h\left\{ {P \cdot \left( {u_{xAcri} - u_{xA}^{C}} \right) \cdot \left\lbrack {1 + \left( {P \cdot \left( {u_{xAcri} - u_{xA}^{C}} \right)} \right)^{\gamma}} \right\rbrack} \right\}} + \frac{1}{2}} = 0.0104}} & {{Eq}\text{.23}} \end{matrix}$

-   -   6. The gradients of the objective function V(ρ) and the         constraint functions g_(xA)(ρ,X,I) and g_(yA)(ρ,X,I) with         respect to the design variables are calculated as follows:     -   6.1) The gradient of the objective function is given by Eq.24:

$\begin{matrix} {\frac{\partial{V(\rho)}}{\partial\rho_{e}} = {\frac{1}{190 \times 140}\left( {{e = 1},2,\cdots,{190 \times 140}} \right)}} & {{Eq}\text{.24}} \end{matrix}$

-   -   6.2) The gradient solution process of the constraint functions         is explained in detail by taking g_(xA)(ρ,X,I) as an example:     -   6.2.1) A gradient expression of g_(xA)(ρ,X,I) is written         according to a chain rule, as shown in Eq.25:

$\begin{matrix} {\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial\rho_{e}} = {{\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{xA}^{L}}\frac{\partial{\overset{\sim}{u}}_{xA}^{L}}{\partial\rho_{e}}} + {\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{xA}^{R}}\frac{\partial{\overset{\sim}{u}}_{xA}^{R}}{\partial\rho_{e}}\left( {{e = 1},2,\cdots,{190 \times 140}} \right)}}} & {{Eq}\text{.25}} \end{matrix}$

-   -   6.2.2) A function in the bracket of tanh(⋅) in Eq.23 is recorded         as R(u_(xAcri)) so the gradient terms ∂g_(xA)(ρ,X,I)/∂ũ_(xA)         ^(L) and ∂g_(xA)(ρ,X,I)/∂ũ_(xA) ^(R) in Eq.25 can be calculated         according to Eq.26:

$\begin{matrix} {\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{xA}^{*}} = {{\frac{{sech}^{2}\left( {R\left( u_{xAcri} \right)} \right)}{2} \cdot \frac{\partial{R\left( u_{xAcri} \right)}}{\partial{\overset{\sim}{u}}_{xA}^{*}}}{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}\text{.26}} \end{matrix}$

-   -    where the gradient term ∂R(u_(xAcri))/∂ũ_(xA) ^(*) is         specifically:

$\begin{matrix} {\frac{\partial{R\left( u_{xAcri} \right)}}{\partial{\overset{\sim}{u}}_{xA}^{L}} = {\frac{u_{xAcri} - u_{xA}^{C}}{2P^{2}} - \frac{P}{2} + \frac{\left( {u_{xAcri} - u_{xA}^{C}} \right)^{\gamma + 1}}{2P^{2}} - {\frac{P\left( {\gamma + 1} \right)}{2}\left( {u_{xAcri} - u_{xA}^{C}} \right)^{\gamma}}}} & {{Eq}.27} \end{matrix}$ $\begin{matrix} {\frac{\partial{R\left( u_{xAcri} \right)}}{\partial{\overset{\sim}{u}}_{xA}^{R}} = {{- \frac{u_{xAcri} - u_{xA}^{C}}{2P^{2}}} - \frac{P}{2} - \frac{\left( {u_{xAcri} - u_{xA}^{C}} \right)^{\gamma + 1}}{2P^{2}} - {\frac{P\left( {\gamma + 1} \right)}{2}\left( {u_{xAcri} - u_{xA}^{C}} \right)^{\gamma}}}} & {{Eq}.28} \end{matrix}$

-   -   6.2.3) The gradient terms ∂ũ_(xA) ^(L)/∂ρ_(e) and ∂ũ_(xA)         ^(R)/∂ρ_(e) in Eq.25 are given in a uniform form according to         the SIMP framework:

$\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{xA}^{*}}{\partial\rho_{e}} = {{- {L_{xA}^{T}\left( K_{xA}^{*} \right)}^{- 1}}\frac{\partial K_{xA}^{*}}{\partial\rho_{e}}{\overset{\sim}{U}}_{xA}^{*}{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}.29} \end{matrix}$

-   -    where K_(xA) ^(*) and Ũ_(xA) ^(*) (*=L,R) are defined in 5.2.2;         and the gradient term ∂K_(xA) ^(*)/∂ρ_(e) is calculated         according to Eq.30:

$\begin{matrix} {\frac{\partial K_{xA}^{*}}{\partial\rho_{e}} = {p\rho_{e}^{p - 1}\left\langle k_{exA}^{*} \right\rangle{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}.30} \end{matrix}$

-   -    where k_(exA) ^(*) is the element stiffness matrix extracted         from K_(xA) ^(*) and         k_(exA) ^(*)         is a square matrix reconstructed by performing a combined         operation on elements in k_(exA) ^(*) according to the element         stiffness matrix, and is consistent with K_(xA) ^(*) in         dimensionality; and     -   6.2.4) All the gradient terms in Eq.26 to Eq.30 are substituted         into Eq.25 to obtain the gradient of the constraint function         g_(xA)(ρ,X,I);

$\begin{matrix} {{\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial\rho_{2}} = 355.29},\ldots,{\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial\rho_{10 \times 70}} = {- 171.04}},\ldots,{\frac{\partial{g_{xA}\left( {\rho,X,I} \right)}}{\partial\rho_{190 \times 140}} = {{- 3}95.63}}} & {{Eq}.31} \end{matrix}$

-   -   7. The design variable is updated by using a moving asymptote         algorithm as follows:

ρ₁=0.76,ρ₂=0.76, . . . ,ρ_(10×70)=0.31, . . . ,ρ_(190×140)=0.78  Eq. 32

The difference between the objective function values in the current and previous iterations is checked, since it is the first iteration, the difference is defined as the current objective function value, which does not meet a convergent threshold value 0.01, and steps 5 to 7 are repeated; and

A finally obtained optimal solution is partially cut out as follows:

ρ₁=1.00,ρ₂=1.00, . . . ,ρ_(10×70)=1E−3, . . . ,ρ_(190×140)=1.00  Eq.33

Iterative optimization converges at the 143^(rd) generation, and the topological structure corresponding to the optimal solution is shown in FIG. 5 ; and the optimal objective performance is V(ρ)=0.4388 while the reliability of constraint performance u_(xA) and u_(yA) at the optimal solution is 0.9920 and 0.9894 respectively, which meet the reliability design requirement for the support structure in the outer cutterhead of the tunnel boring mechanism. After further smoothing the contour of the result of reliability-based topology optimization design, the final design of the support structure in the outer cutterhead of the tunnel boring mechanism is shown in FIG. 6 .

It is to be noted that the content and specific implementations of the present disclosure are intended to prove the practical application of the technical solution provided by the present disclosure, and should not be interpreted as limiting the scope of protection of the present disclosure. Any modification and change of the present disclosure within the scope of protection of the spirit of the present disclosure and the claims fall within the scope of protection of the present disclosure. 

What is claimed is:
 1. A reliability-based topology optimization design method for a part structure considering bounded hybrid uncertainties, comprising: step 1) considering following uncertainties in manufacture and service processes of the part structure: regarding an amplitude and a direction of an external load with insufficient sample information as interval uncertainties, regarding a material property of the part structure with sufficient sample information as a bounded probabilistic uncertainty, and describing a bounded probabilistic uncertainty parameter as a random variable subject to generalized beta distribution; step 2) discretizing a design domain of the part structure, comprising: simplifying a force condition of the part structure into a two-dimensional plane stress state, retaining installation holes and removing structural details to improve calculating efficiency, placing the simplified part structure in a regular rectangular design domain, dividing the rectangular design domain into N_(x)×N_(y) square elements, where N_(x) and N_(y) are numbers of divisions along x, y axes, respectively; and imparting, based on a solid isotropic material with penalization (SIMP) topology optimization framework, each element with a unique design variable ρ_(e)∈[0,1] (e=1, 2, . . . , N_(x)·N_(y)); step 3) imposing physical constraints and geometric constraints on the discretized structure, comprising: step 3.1) imposing the physical constraints including fixing or supporting and an external load according to a classical finite element method; and step 3.2) the geometric constraints including specified holes in the structure and areas where materials are forcibly retained, setting design variables corresponding to elements in the holes as ρ_(e)≡0, setting design variables corresponding to elements in the areas where the materials are forcibly retained as ρ_(e)≡1, and keeping values of the design variables corresponding to elements in the holes and the areas where the materials are to be retained unchanged in subsequent optimization process; step 4) establishing, by taking a space utilization rate of the design domain as an objective function and taking displacements of several key points of the part structure under a joint influence of interval and bounded probabilistic hybrid uncertainties as reliability constraint performance, a reliability-based topology optimization design model for the part structure, as shown in Eq.1: $\begin{matrix} {{\min\limits_{\rho}{V(\rho)}} = {\sum\limits_{e = 1}^{N_{e}}{\rho_{e}/V_{0}}}} & {{Eq}.1} \end{matrix}$ s.t.g_(q)(ρ, X, I) = −P(u_(q)(ρ, X, I) ≤ u_(qcri)) + P_(q) ≤ 0 (q = 1, 2, …, N_(con)) K(ρ, X)U = F(I) 0 ≤ ρ_(min) ≤ ρ_(e) ≤ 1(e = 1, 2, …, N_(e)) X = (X₁, X₂, …, X_(m))^(T), I = (f₁, f₂, …, f_(n), α₁, α₂, …, α_(n))^(T) where in Eq.1, ρ=(ρ₁, ρ₂, . . . , ρ_(N) _(x) _(·N) _(y) )^(T) is a design vector composed of design variables ρ_(e)(e=1, 2, . . . , N_(x)·N_(y)), and ρ_(min) is a minimum allowable value for the design variables; a total number of elements is N_(e)=N_(x)·N_(y); a bounded probabilistic uncertainty vector X=(X₁, X₂, . . . , X_(m))^(T) comprises m uncertain material properties of the part structure; an interval uncertainty vector I=(f₁, f₂, . . . , f_(n), α₁, α₁, . . . , α_(n))^(T) comprises amplitudes f₁, f₂, . . . , f_(n) and direction angles α₁, α₂, . . . , α_(n) of n uncertain external loads on the part structure; where V(ρ) is the space utilization rate of the design domain, corresponding to a total material usage of the part structure, and V₀ is a volume of the design domain; where g_(q)(ρ,X,I) is a q^(th) constraint function; u_(q)(ρ,X,I) is a displacement at a q^(th) key point regarded as a q^(th) constraint performance, denoted as u_(q) for short, and u_(qcri) is an allowable maximum value of u_(q)(ρ,X,I); P(⋅) calculates an occurrence probability of an event in a bracket, P_(q) is a reliability index of the q^(th) constraint performance, and N_(con) is a number of constraint functions; and where in an equilibrium equation K(ρ,X)U=F(I) of the part structure, K(ρ,X) is a (2(N_(x)+1)(N_(y)+1))×(2(N_(x)+1)(N_(y)+1))-dimensional global stiffness matrix, which is affected by the design vector ρ and the bounded probabilistic uncertainty vector X; F(I) is a (2(N_(x)+1)(N_(y)+1))-dimensional nodal force vector affected by the interval uncertainty vector I; U is a (2(N_(x)+1)(N_(y)+1))-dimensional nodal displacement vector; and u_(q) is extracted from U according to u_(q)=L_(q) ^(T)U, wherein in a (2(N_(x)+1)(N_(y)+1))-dimensional column vector L_(q), the element at a position corresponding to the q^(th) key point is 1, and other elements are all 0; step 5) calculating reliability of a constraint performance under bounded hybrid uncertainties, comprising: step 5.1) searching a worst working condition of constraint performance u_(q), comprising: step 5.1.1) setting X=μ_(X)=(μ_(X) ₁ , μ_(X) ₂ , . . . , μ_(X) _(m) ) where μ_(X) ₁ , μ_(X) ₂ , . . . , μ_(X) _(m) are means of uncertainties X₁, X₂, . . . , X_(m), the constraint performance u_(q) is only affected by the interval uncertainty vector; and an uncertain external load F_(s) (s=1, 2, . . . , n) is rewritten as its components in horizontal and vertical directions such as F_(s)=[f_(s) cos α_(s),f_(s) sin α_(s)]^(T); step 5.1.2) calculating, based on a linear elastic hypothesis in small deformation, displacement U_(s) caused by F_(s) according to Eq.2 based on e_(s) ^(x)=[1 0]^(T), e_(s) ^(y)=[0 1]^(T) and F_(s), where e_(s) ^(x) and e_(s) ^(y) are respectively element nodal forces in the horizontal and vertical directions at an point where F_(s) exerts; U _(s) =U _(s) ^(x) +U _(s) ^(y) =u _(s) ^(x) f _(s) cos α_(s) +u _(s) ^(y) f _(s) sin α_(s) =[u _(s) ^(x) ,u _(s) ^(y) ]·F _(s)  Eq.2 where u_(s) ^(x)=[u_(s) ^(x) 0]^(T) and u_(s) ^(y)=[0 u_(s) ^(y)]^(T) are respectively nodal displacement vectors calculated through the equilibrium equation of the part structure when only the element nodal force e_(s) ^(x) or e_(s) ^(y) acts; step 5.1.3) calculating gradients of the constraint performance u_(q) with respect to an amplitude and a direction of an uncertain external load according to Eq.3 and Eq.4: $\begin{matrix} {\frac{\partial u_{q}}{\partial f_{s}} = {{L_{q}^{T}\frac{\partial U}{\partial f_{s}}} = {{L_{q}^{T}{\sum\limits_{i = 1}^{n}\frac{\partial U_{i}}{\partial f_{s}}}} = {L_{q}^{T} \cdot \left\lbrack {u_{s}^{x},u_{s}^{y}} \right\rbrack \cdot \left\lbrack {{\cos\alpha_{s}},{\sin\alpha_{s}}} \right\rbrack^{T}}}}} & {{Eq}.3} \end{matrix}$ $\begin{matrix} {{\frac{\partial u_{q}}{\partial\alpha_{s}} = {{L_{q}^{T}\frac{\partial U}{\partial\alpha_{s}}} = {{L_{q}^{T}{\sum\limits_{i = 1}^{n}\frac{\partial U_{i}}{\partial\alpha_{s}}}} = {L_{q}^{T} \cdot \left\lbrack {u_{s}^{x},u_{s}^{y}} \right\rbrack \cdot \left\lbrack {{{- f_{s}}\sin\alpha_{s}},{f_{s}\cos\alpha_{s}}} \right\rbrack^{T}}}}};} & {{Eq}.4} \end{matrix}$  and step 5.1.4) solving, by utilizing results of Eq.3 and Eq.4, a worst working condition Ĩ_(q) by a gradient search algorithm, and an external load of the worst working condition being F={tilde over (F)}_(q); step 5.2) restoring μ_(x) to X under the worst working condition Ĩ_(q), so that constraint performance ũ_(q)=u_(q)(ρ,X,Ĩ_(q)) under the worst working condition is only manifested as a probabilistic type, and solving a performance fluctuation under the worst working condition to evaluate reliability, comprising: step 5.2.1) calculating, according to Eq.5, a gradient of the constraint performance ũ_(q) with respect to a bounded probabilistic uncertainty parameter X_(i)(i=1, 2, . . . , m) under the worst working condition: $\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{q}}{\partial X_{i}} = {{L_{q}^{T}\frac{\partial\overset{\sim}{U}}{\partial X_{i}}} = {{{- L_{q}^{T}}K^{- 1}\frac{\partial K}{\partial X_{i}}\overset{\sim}{U}} = {{- L_{q}^{T}}{K^{- 1}\left( {\sum\limits_{e = 1}^{N_{e}}{\rho_{e}^{p}\frac{\partial k_{e}}{\partial X_{i}}}} \right)}\overset{\sim}{U}}}}} & {{Eq}.5} \end{matrix}$ where a summation symbol is a combination operation of the element stiffness matrix defined by a finite element theory, k_(e) is an element stiffness matrix, a nodal displacement vector Ũ under the worst working condition is obtained by solving a governing equation KŨ={tilde over (F)}_(q) under the worst working condition; and p is a penalty factor; step 5.2.2) searching, based on a result of Eq.5, two bounded probabilistic uncertainty vectors minimizing or maximizing u_(q)(ρ,X,Ĩ_(q)) respectively according to Eq.6: $\begin{matrix} \left\{ {\begin{matrix} {X_{q}^{L} = {\underset{X}{argmin}{u_{q}\left( {\rho,X,{\overset{\sim}{I}}_{q}} \right)}}} \\ {X_{q}^{R} = {\underset{X}{argmax}{u_{q}\left( {\rho,X,{\overset{\sim}{I}}_{q}} \right)}}} \end{matrix}\left( {{q = 1},2,\ldots,N_{con}} \right)} \right. & {{Eq}.6} \end{matrix}$ where corresponding to X_(q) ^(L) and X_(q) ^(R), global stiffness matrix K is K_(q) ^(L) and K_(q) ^(R), respectively, K_(q) ^(L) and K_(q) ^(R) are uniformly marked as K_(q) ^(*)(*=L,R), the nodal displacement vector Ũ is Ũ_(q) ^(L) and Ũ_(q) ^(R), respectively, and Ũ_(q) ^(L) and Ũ_(q) ^(R) are uniformly marked as Ũ_(q) ^(*)(*=L,R); and step 5.2.3) denoting ũ_(q) ^(L)=u_(q)(ρ,X_(q) ^(L),Ĩ_(q)) and ũ_(q) ^(R)=u_(q)(ρ,X_(q) ^(R),Ĩ_(q)), defining [ũ_(q) ^(L),ũ_(q) ^(R)] as a performance fluctuation of the constraint performance u_(q) under the worst working condition, and calculating a reliability {tilde over (R)}_(q) of the constraint performance according to Eq.7: $\begin{matrix} {{\overset{\sim}{R}}_{q} = {{\frac{1}{2}\tanh\left\{ {P \cdot \left( {u_{qcri} - u_{q}^{C}} \right) \cdot \left\lbrack {1 + \left( {P \cdot \left( {u_{qcri} - u_{q}^{C}} \right)} \right)^{\gamma}} \right\rbrack} \right\}} + \frac{1}{2}}} & {{Eq}.7} \end{matrix}$ where ũ_(q) ^(C)=(ũ_(q) ^(R)+ũ_(q) ^(L))/2 is a midpoint of the performance fluctuation under the worst working condition; a multiplier is P=1/(u_(q) ^(W)−ε^(u)) where ε^(u) is a small constant for adjusting reliability at a boundary position of the performance fluctuation under the worst working condition; u_(q) ^(W)=(ũ_(q) ^(R)−ũ_(q) ^(L))/2 is a radius of the performance fluctuation under the worst working condition; and γ∈{2i|i∈N⁺} is a regulatory factor; and establishing the constraint functions in Eq.1 as Eq.8: g _(q)(ρ,X,I)=−{tilde over (R)} _(q) +P _(q)≤0 (q=1,2, . . . ,N _(con))  Eq.8 step 6) calculating gradients of objective and constraint functions with respect to the design variables: step 6.1) calculating the gradient of the objective function through Eq.9: $\begin{matrix} {\frac{\partial{V(\rho)}}{\partial\rho_{e}} = {\frac{1}{V_{0}}\left( {{e = 1},2,\ldots,N_{e}} \right)}} & {{Eq}.9} \end{matrix}$ step 6.2) solving the gradient of a constraint function as follows: step 6.2.1) writing a gradient expression of g_(q)(ρ,X,I) according to a chain rule, as shown in Eq.10: $\begin{matrix} {\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial\rho_{e}} = {{\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{q}^{L}}\frac{\partial{\overset{\sim}{u}}_{q}^{L}}{\partial\rho_{e}}} + {\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{q}^{R}}\frac{\partial{\overset{\sim}{u}}_{q}^{R}}{\partial\rho_{e}}\left( {{e = 1},2,\ldots,\ N_{e}} \right)}}} & {{Eq}.10} \end{matrix}$ step 6.2.2) denoting a function in a bracket of tanh(⋅) in Eq.7 as R(u_(qcri)), and calculating gradient terms ∂g_(q)(ρ,X,I)/∂ũ_(q) ^(L) and ∂g_(q)(ρ,X,I)/∂ũ_(q) ^(R) in Eq.10 according to Eq.11: $\begin{matrix} {\frac{\partial{g_{q}\left( {\rho,X,I} \right)}}{\partial{\overset{\sim}{u}}_{q}} = {\frac{\sec{h^{2}\left( {R\left( u_{qcri} \right)} \right)}}{2} \cdot \frac{\partial{R\left( u_{qcri} \right)}}{\partial{\overset{\sim}{u}}_{q}^{*}}{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}.11} \end{matrix}$ where a gradient term ∂R(u_(qcri))/∂ũ_(q) ^(*) is as follows: $\begin{matrix} {\frac{\partial{R\left( u_{qcri} \right)}}{\partial{\overset{\sim}{u}}_{q}^{L}} = {\frac{u_{qcri} - u_{q}^{C}}{2P^{2}} - \frac{P}{2} + \frac{\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma + 1}}{2P^{2}} - {\frac{P\left( {\gamma + 1} \right)}{2}\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma}}}} & {{Eq}.12} \end{matrix}$ $\begin{matrix} {\frac{\partial{R\left( u_{qcri} \right)}}{\partial{\overset{\sim}{u}}_{q}^{R}} = {{- \frac{u_{qcri} - u_{q}^{C}}{2P^{2}}} - \frac{P}{2} - \frac{\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma + 1}}{2P^{2}} - {\frac{P\left( {\gamma + 1} \right)}{2}\left( {u_{qcri} - u_{q}^{C}} \right)^{\gamma}}}} & {{Eq}.13} \end{matrix}$ step 6.2.3) giving gradient terms ∂ũ_(q) ^(L)/∂ρ_(e) and ∂ũ_(q) ^(R)/∂ρ_(e) in Eq.10 in a uniform form according to the SIMP framework: $\begin{matrix} {\frac{\partial{\overset{\sim}{u}}_{q}^{*}}{\partial\rho_{e}} = {{- {L_{q}^{T}\left( K_{q}^{*} \right)}^{- 1}}\frac{\partial K_{q}^{*}}{\partial\rho_{e}}{\overset{\sim}{U}}_{q}^{*}{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}.14} \end{matrix}$ where K_(q) ^(*) and Ũ_(q) ^(*) (*=L,R) are defined in 5.2.2; and a gradient term ∂K_(q) ^(*)/∂ρ_(e) is calculated according to Eq.15: $\begin{matrix} {\frac{\partial K_{q}^{*}}{\partial\rho_{e}} = {p\rho_{e}^{p - 1}\left\langle k_{eq}^{*} \right\rangle{\text{(*}\left. {{= L},R} \right)}}} & {{Eq}.15} \end{matrix}$ where k_(eq) ^(*) is an element stiffness matrix extracted from K_(q) ^(*), and

k_(eq) ^(*)

is a square matrix reconstructed by performing a combined operation on elements in k_(eq) ^(*) according to the element stiffness matrix, and is consistent with K_(q) ^(*) in dimensionality; and step 6.2.4) substituting all the gradient terms in Eq.11 to Eq.15 into Eq.10 to obtain a gradient of the constraint function g_(q)(ρ,X,I); and step 7) updating the design variables by using a moving asymptote algorithm based on the gradients of the objective and constraint functions with respect to the design variables, checking a difference value between an objective function value in a current iteration and an objective function value in a previous iteration, wherein the difference value for a first iteration is defined as an objective function value, and when the difference value is less than a convergence threshold, outputting the updated design variables, and otherwise, repeating the steps 5) to 7). 